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1, INTRODUCTION 

Nowadays, PID controllers have received considerable attention in the last years both from an 
academic and industrial point of view [1-5]. In fact, in principle, they provide more flexibility in the 
controller design, concerning the standard PID controllers, because they have five parameters to select. 
However, this also implies that the tuning of the controller can be much more complicated. They have been 
successfully applied in practical applications such as motion control of manipulators and chaos control of 
electrical circuits. In these applications, it has been verified that PID controllers can improve the performance 
of traditional control system adopting integer order PID controllers. The most important advantage of the PID 
controllers is that they can afford more extensive possibilities offered by their additional fractional order 
dynamics [4-7]. However, this also indicates that the tuning strategies of PID controllers are much more 
complicated. In the researches on the PID controllers, tuning of controller parameters has become a 
significant issue. 

In general, the tuning methods for FID controllers are classified into analytical, numerical, and rule- 
based ones. In [6-7] the controller parameters have been analytically derived by solving nonlinear equations 
fulfilling the gain/phase crossover frequency and phase/gain margin specifications. The robustness to loop 
gain variations specification proposed in [8] has also been widely used to design FID and proportional— 
integral (PI) controllers. The merits of the analytical method are obvious; however, it is available only when 
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the equations are few in number and simple. Therefore, it is very difficult to obtain a complete FID controller 

for the MB system by solving five complicated nonlinear equations. As for the rule-based method, it can 

easily calculate the controller parameters based on empirical tuning rules, which can be observed in [9-12]. 

In this paper, we propose and investigate optimal tuning PID controller of unstable fractional order system by 

desired transient characteristics using real interpolation method (RIM). The main advantages of this method 

are drawn as the followings: 

a) This is the novel practical method based on the desired settling time and overshoot percentage. 

b) The results are close to the desired parameters. 

c) The novel method can tune an unstable fractional order system by real interpolation method (RIM). 

d) The novel method is simplicity and computer efficiency. 

e) The novel method can find an optimal solution for tuning task in both academic and industrial purposes. 
The paper is organised as follows. In Section 2 the problem is formulated. The tuning rules are 

described in Section 3. Simulation results are presented in Section 4, where a comparison with other tuning 

rules is performed. Conclusions are drawn in Section 5. 


2. RESEARCH METHOD 
2.1. Real Interpolation Method (RIM) 

The real interpolation method is one of the methods, which works on mathematical descriptions of 
the imaginary area. The method is based on real integral transform as follow [13]: 


F(8)= [ f(t)e~Stdt,5 €[C,00),C 20 
0 (1) 


which assigns to the original f(t) the image F(6) as a function of the real variable ° Formula of direct 
transform (12) can be considered as a special case of the direct Laplace transform by replacing the complex 
variable p = 6 + jw for real 6 variable. Another step towards the development of the instrumentation 
method is the transition from continuous functions F(6) to their discrete form, using the computing resources 
and numerical methods. For these purposes, real interpolation method is represented by numerical 
characteristics {F (6; )},. They are obtained as a set of values of the function F(6) in the nodesi = 1,7, 
where 77 is the number of elements numerical characteristics, called its dimension. 

Selecting of interpolations 6; is a primary step in the transition to a discrete form, which has a 
significant impact on the numerical computing and accuracy of problem solutions. Distribution of nodes in 
the simplest variant is uniform. Another important advantage of real interpolation method is cross-conversion 
properties [11]. It dues to the fact that the behavior of the function F(6) for large values of the argument 6 is 
determined mainly by the behavior of the original f(t) for small values of the variable t. In the opposite 
case, the result is the same: the behavior of the function F(é6) for small values of the argument 6 is 
determined mainly by the behavior of the original f(t) for large values of the variable t [13]. 

When considering the original f(t) of dynamic characteristics of dynamic systems, formula (1) 
leads to an operator model, which under certain conditions can be considered as special cases of the models 
based on the Laplace transform. Thus, in (1) replacing of the function f(t) by k(t) - the impulse transient 
function of the dynamic system, we obtain its transfer function. From here we can find the elements of a 
discrete model of the system, and its transfer function by performing the discretization procedures for nodes 
6;,t€ 1,7: 








W(5)= fe(t)e Mar, teh 
0 (2) 


Function W(6) is a real transfer function of control automatic systems, having an impulse transient 
response k(t). Function W(6) could be received based on determination of transfer function such as a 
relationship of the imaginary of output Y(6) and X(6) input signals 


W(6)=Y(6)/X(6) (3) 


In which the imaginary of the output signal and the input signal is calculated from the original 
functions of the input x(t) and output y(t) signals: 
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[o.0) 


Y(6) =|(e)e at 


(4) 


X (6) =|x(t)e "adr 
: (5) 
where x(t) input signal, and y(t) output signal of the system. 
The input-output relationship of the system has a form like below: 
¥ (5) =W(5)X(8) S 
2.2. Problem formulation 
In this section, the control system with negative unity feedback is drawn in Figure 1. 


u(t) y(t) 








Controller 


Figure 1. Control system with negative unity feedback 


Transfer function of PID controller can be formulated as. 


C(s) = ca K,s 
§ (7) 


Transfer function of plant (in fractional form) 1s 


GQ B,(s) _ bs +b sh 4--4+bs" +b s” 

A,(s) 4,5°" +4,_,8°™! +---+a,57 +a,5” (8) 
where m,n € N,and Qo, ...,; Am, Do) ++, Dn, Am > Am_-1 > > Ay > A = 0,8, > Pa_1 > > > Po = 
0 are arbitrary real numbers a,, # 0,b,, # 0. 

Transfer function of PID controller (7) can be rewritten in the rational form as the following 


B-(s)_ K,s°+K,s+K, 


ee A. (Ss) KY 





(9) 


Characteristic equation 1s 


W(s) =1+ C(s)G(s) (10) 


Characteristic polynomial is 


P(s) = A.(s)Ag(s) + B.(8) Bg (8) (11) 


A gain-phase margin tester (GPMT) can be thought of as a “virtual compensator”, provides 
information for plotting the boundaries of constant gain margin and phase margin in a parameter plane. 
The frequency independent GPMT is given in the form [15]: 


G,(M,¢) = Me”. (12) 
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For a given IOPID controller parameters os es a the closed-loop system is said to be bounded- 


P(s, K,,,K;,,K,) 


input bounded-output (BIBO) stable if the quasi-polynomial has no roots in the closed right- 


half of the s-plane (RHP). The stability domain S in the parameter space P with ae ae being coordinates 


is the region that for oo the roots of quasi-polynomial PoE ta) all lie in open left-half of 


the s-plane (LHP). The boundaries of the stability domain S which are described by real root boundary 
(RRB), infinite root boundary (IRB) and complex root boundary (CRB) can be determined by the D- 


P(O,K)=0 P(o,K)=0 


decomposition method [14]. These boundaries are defined by the equations and 


P(EJ@.K)=0 fo-@€(0,%) respectively, where PCS) 


system and K the vector of controller parameters. 
a) Determining RRB 

In applying the descriptions of stability boundaries of the stability domain S to the FOCE in (6), 
the RRB turns out to be simply a straight line given by 


is the characteristic function of the closed-loop 


P(0,K,,K,,K,) =0=> K, =0, (13) 


for 8” =1 jn the transfer function of the plant in (8). 
b) Determining IRB 

There is more theoretical difficulties for the calculating of the IRB due to fractional component. 
FOCE possesses an infinite number of roots, which cannot be calculated analytically in the general case. 
However, the asymptotic location of roots far from the origin is well known [21], [22], which may lead to 
IRB. The objective of this section is to determine the stabilizing region in (Kp, Ki) plane with given Kd, 
and values for which the following complex polynomial is 


D(s) = sL(s)A,(s)+(K,s° +K,,s+K,)M(s)B,(s) (14) 


where L(s) and M(s) are given complex fractional order polynomials. When L(s) = M(s) = 1, the stabilization 
of (14) reduces to the standard IOPID stabilization. 
CO ea Suppose when *~ % , [L(s)Ao(s) ]/[M (5) Be (s) Js where t and c are real and complex 
numbers respectively. We have: 
(a) If 7~!, then the boundary does not exist. 
(b) If =! andc isnot real, then the boundary does not exist. 
(c) If 7=! andc is real, then 
N(s)(Kps+Kj+Kqs*) 


D(s) __ iy. 
pee sL(s)Ag(s) — lim 1+ sSAgc(s) 


=14+K,c=0 (15) 


Ko" (16) 


Cc 


c) Determining CRB 
To construct the CRB, we substitute *~/@ into (11) to obtain, 


PG) = A.Ga@)A, Go) + B.Ga@)B,G@) =0 (17) 


Using D-decomposition to find stability [21], to construct the CRB, we substitute *~ J© into (11) to obtain 


Using D-decomposition to find stability [ [1] Reference -A graphical tuning] 
j? =e?” =cos(~g) + jsin(¢) 
In s plane 2 2, where “is a real number. The formula (12) can be 
expressed into a formula with the separate real and imaginary components 
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PQG@) = | (K,- K, O° )R pg (@) — Oly (@) | * jo| K Rac (@) + R ,<(@) | =O (16) 


where Roc Rig the are the real component of Pao) the real and imaginary components of AcQ@) 
respectively. Finally, by setting the real and imaginary parts equalized zero the formula (15) leads to the 


detail form: 


oe +K R,¢(@) =0 


(K,—-K, @)R,,(@) — @1,,(@) =9 (17) 


The CRB is determined by formula (17). 


2.3. Algorithm 
2.3.1 Design of the desired transient characteristics 

Performance qualities of the control system can be evaluated by settling time, overshot percentage, 
damping ratio, natural frequency, rise time. There are some methods to design the desired transfer 
characteristics [18-20]. In this paper, the desired transient characteristics could be used a second order 
system. A second order system can be described by the following form without zeros 


2 


a) 
H(s) = ————_—_; 
s°+2o0 s+; (18) 
With zeros, 
2 
H(s)= so, /(ag)+ a; 


2 


s°+260,5+ 0° (19) 


Where “ natural frequency, S damping ratio, % coefficient related to overshoot percentage. And there is a 
relationship between settling time and natural frequency and damping ratio 


be tS 


cat 
set O, é (20) 


Another option is that we can use the mat lab or another software program to design a desired 
transient characteristics. 


2.3.2 Real interpolation method 
The desired transient characteristics could be given by massive data, table and transient function. 


Let see the transient characteristics given by form: function MmO=FO or massive "0 =F 
From Figure 1, the equivalent transfer function of the feedback system has a form 


W(5) = LOG) 
1+C(s)G(s) (21) 


Or it could be rewritten by input and output transfer function 


Ys) __ C(S)G() 
X(s) 1+C(s)G(s) (22) 


X (Ss) 


If the input of the tuning system is a step function, when is determined 


X(s)=1/s (23) 
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To determine coefficient Kp, Ki, Kd in PID controller for fractional order transfer function using the 
real interpolation method. By substituting *=° formula (22) can be rewritten in the real transfer form. 


Y(0) 2 C(0)G(6) 
X(5)_ 1+C()G(S) (24) 


From (24) we obtain the formula determined PID controller. 


1 


( ———————— 
G(6)[X(6)/ Y(6) -1] (25) 


where [(O) X(O) Y(O) 
respectively. 


real transfer function of the plant, the imaginary of output and input signals 


X(0), ¥(O) could be determined by formula (4) and (5) in the real form or by 


X(0), Y(O) 


The functions 
following formula. If the transient characteristics are given by a massive the functions can be 


obtained by the following formula: 


x (6, ) = x(t,)e°" At 
i=l (26) 


ZB (27) 


t. h(t, 
where «(1) and (4) are the input value and the output value at time 7 respectively; 4? -the 
sample time and N - the number of samplings. 


X(0) 


If the input of the tuning system is a simple step function, when is determined in real form. 


X(6)=1/6 (28) 
The number of the unknown coefficients of the PID controller ( ”° ' ) is three; it means that 


the dimension of the elements numerical characteristics is three (7 ~ ).Selecting value of nodes is a primary 


t (s 
step. It effects on the accuracy of the tuning task. The meaningful region is from to Se ~~ . It means that 
the maximum value of first node % can be defined by the following condition: the value of the integral in (4) 
by the settling time ty reduces to a negligible value 4 = 9-001+0.05 which satisfy the condition 
Ol ser < : 
h(t, )e™ <A . Hence calculated expression for the node ©. can be shown below. 


_In(A/h(t,,,)) 


5,= 
Eset (29) 


where A -calculation error, ME ser) settling value, ‘set settling time. 
Other node in the simplest case, can be determined by the nominal distribution. 


0, =10,,1=2,n (30) 


In case, the result according to the maximum value of the first node do not meet the desired output 
performances, the value of the first node should be decreased. In this paper, by varying the value of the first 
node from the maximum, determined by (29) approach to 0, we carry out the optimization of the range of 
nodes for the output requirements of the desired system. 

The numerical characteristics of the PID controller’s real transfer function. 
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C(6,)=K,,+ . +K,6, i=1,7 








(31) 
Oe, 1 
; aos K_ ,K.,K 
We have a system of equations with unknown coefficient ” “ ¢ 
115 
O, 
eal: 28 
5, K, F(6,) 
K =| K, F =| Fé 
1 a3 b, 1 ( 5) 
O; K, F(0;) (33) 
ae K KK 
The unknown coefficients ”’ ‘’ “% can be obtained in matrix form 
K=C'F (34) 
Algorithm: 
a) Determination of stability region (13), (16) and (17). 
b) Design a desired transient characteristic (18) or (19) 
c) Calculation value of nodes (29) and numerical characteristics (22) or (28), (27), (31) and (32). 
d) Solving the equation to find the parameters of PID controller (33) and (34). 
e) Estimation of the tuned PID controller in the time and frequency domains. 
f) Investigation of the optimal range of the values of node by varying the value of the nodes. 
3. NUMERICAL RESULTS AND DISCUSSION 
Given fractional order transfer function of a unstable bearing system [zhong16], 
oe 6438. 
s*7° + 330.045'*! — 84268.855°** —15869154.14 (35) 


Investigation of STABILITY REGION of the unstable fractional order transfer function of the 
bearing system (11). 
a) Determine RRB 

It 1s important to determine the stability region of the PID controller, tuned for the system (35). 


The RRB can be determined by (13) as the bellow 


b) Determine IRB 
The IRB can be determined according to (16) 





en 2OUS) te ee es, 
msD(s) ) 
Not exist the IRB. 


c) Determine CRB 

The CRB can be determined by three conditions (18), (19) and (20). The results are shown in the 
following figures. The following figure shows the stability region of the transfer function (11). The CRB in 
(Kp, Ki) plane when Kd=10 is shown in the Figure 2. The cross-line area is a feasible region for the 
coefficients of the PID controller. 
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Ki 








7 | | 0 a ole 6000 oh oly 9000 10000 
Figure 2. The feasible region of coefficients Kp, Ki, Kd 
The cross-line region demonstrates stability region of Kp, Ki where Kd=10. Oviously, the stability 


region of PID controller’s coefficients lie in the first quarter. In the Figure 3 there are shown more feasible 
area with five values of Kd. 
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Figure 3. Stability region of PID controller depending on its coefficients 


All boundaries (RRB, IRB and CRB) are shown in following figure. Kp, Ki, Kd stability region 
The results simply lead to the conclusion 


IK >= 2465 


K. 
<0 (39) 
K,,K,,K, ; 
Three parameters of ” ' are located in the first quarter of coordinate axis system. 
STUDY CASE 1 
Design of the desired transient characteristics. 
Using (22) and (23) with the desired transient time Tr =9-03 8 and overshoot POT = 50% we 


formed a transfer function with the step response shown in the Figure 4. According to Figure 4, the transient 


time of the desired system accounts for Te = 0.0306 8 and the overshoot is given a desired value at 


POT = 50.1% 
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System: Wd 

Peak amplitude: 1.5 
Overshoot (%): 50.1 
At time (seconds): 0.00937 


step Response 








Amplitude 





O 0.01 0.02 0.03 0.04 0.05 0.06 
Time (seconds) 


Figure 4. The desired transient characteristics with overshoot percentage 50.1% and settling time 0.0366 s 


Obviously, the number of the numerical characteristics is three (7 ~ 3 ) because of three unknown 


coefficients of the PID controller. Calculation of the nodes of numerical characteristics: 

Using the formula (30) with the input parameters ‘= 0.056 s, At=O0.001s, A=0.001, he = the first 
node can be determined ® =!88-73® Other nodes are determined by the condition of uniform distribution: 
6, = 16, i 5, =[188.736 377.473 566.209] 


1 





sae Finally, the value of nodes can show the following massive 


— F. =}-25230 402.219 6006 
Based on the finding nodes the value of massive F can calculated and make up ' | 2 | 


and C massive accounts for 


1 5.298-10° 188.736 39271.902 
C=|1 2.649-10° 377.473 K =| -1.134-10' 
1 1.7660-10° 566.209 -23.375 


and the coefficients of the PID controller are 


According to the feasible area (39) of PID controller the values of the tuned parameters of the 
controller are not satisfied. When d1={2.5, 5, 5, 7.5, 10 the results are shown in the Figure 5. 





X: 0.009366 Step Response 

















Amplitude 




















1 
0 0.01 0.02 0.03 0.04 0.05, 0.06 0.07 0.08 0.09 
Time (seconds) 





Figure 5. Time response of the tuning systems 
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In terms of the percentage overshot, the tuning systems are close to the desired system. With 61=2.5 
the percentage overshot of the tuning system is slightly higher than the expected value. With other values of 
ol the figure is lower than that of the expected value with the higher value of the first node. According to the 
settling time, the values of the tuning systems are greater than that of the desired system. 

On the other hand, the shape of the time response figures do not following the desired transient 
response. The bode graph show in the following figure. It 1s clear that the bode graphs of the tuning systems 
are close to each other. The magnitude and phase margin of the tuning system are close to each other. 
The bode diagram of the open loop system is shown in Figure 6. 


Bode Diagram 











Magnitude (dB) 











775 
eee 








Phase (deg) 











Ps TTT TTT TT TTT 
10° 10 10- 10 10 Frequency ofrad/s) 10 10° 10° 10° 10° 





Figure 6. Bode of the open loop system 


The detail values of the tuning system are presented in the following table. 
STUDY CASE 2 

We carry out the investigation for the higher value of the desired overshoot percentage and settling 
time. The desired overshoot percentage and settling time account for 100% and 0.0429 s respectively. The 
desired transient process is shown in the Figure 8. The desired transient characteristics with overshoot 
percentage 99.1% and settling time 0.0429 s. 

The results of the tuning process will illustrate in the Figure 7. When 61=5, 10, 20, 50 the time 
responses are shown in the following Figure 7. 





System: Gd Step Response 
Time (2econds): 0.006852 T 
Aumplitude: 1 Seo 





Fl leet ag a i SSIS Sey Seni Sie Sits Sie Sr ei ee eis aia ats 
- 600 s + 27889 
a a hag a pS 
I s* + 250.5 $+ 27889 
ee eee ita! asddemie aeieiescics, “aca esac alate aa. Wiss ace as enc ie acct | ccs elec 
a 
a oe ' 
a 1! ‘ 
= log System: Gd 
=> » ot Settling time (=econds): 0.0425 
= TES CSR RC RCS ISSR OS eee SS SS SS 
=v uae 
I I 
I I 
I I 
ty ee eens a pean ee aaa anise e ab bee Boe eel ater eae a ea terete erg Ses 
I I 
I I 
I I 
0 
0 0.01 0.02 0.03 0.04 0.05 0.06 


Time (seconds) 


Figure 7. The desired transient characteristics with overshoot percentage 99.1% and settling time 0.0429 s 


Indonesian J Elec Eng & Comp Sci, Vol. 14, No. 3, June 2019: 1177-1188 


Indonesian J Elec Eng & Comp Sci ISSN: 2502-4752 O 1187 


Step Response 























Amplitude 

















; 
0 0.05 Time (seconds) ea a 


Figure 8. Time response of the tuning systems 


The shapes of the transient response of the tuning systems are closer to the desired system than that 
in the first example. Higher accuracy than that of the first example. Bode graph of the tuning systems is 
shown in the Figure 9. 
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Figure 9. Bode graph of the tuning systems 


The magnitude margins show the same value for all tuning system, however, the phase magnitude 
margins account for higher value with the higher value of the first node. 


4. CONCLUSION 

In this paper, we propose the method of tuning conventional PID controller for unstable transient 
characteristics. The results show that: 1) This is the novel practical method based on the desired settling time 
and overshoot percentage; 2) The results are close to the desired parameters; 3) The novel method can tune 
an unstable fractional order system by real interpolation method (RIM); 4) The novel method 1s simplicity 
and computer efficiency; 5) The novel method can find an optimal solution for tuning task in both academic 
and industrial purposes. 
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